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Abstract:  Scattering  of  electromagnetic  fields  from  ob- 
jects on  or  beneath  the  ground  can  be  analyzed  with  the 
Finite- Difference  Time-Domain  (FDTD)  method.  The 
equivalence  principle  can  be  utilized  for  generation  of  in- 
cident plane  waves  where  the  Fresnel  reflection  and  re- 
fraction coefficients  are  multiplied  by  the  pulse  spectrum 
in  the  frequency  domain  and  the  result  is  transformed 
into  the  time-domain.  Errors  will  appear  in  the  Huy- 
gens} sources  due  to  numerical  dispersion  where  the  main 
errors  come  from  the  use  of  analytical  reflection  and  re- 
fraction coefficients . Their  influence  on  the  accuracy  is 
generally  negligible  for  applications  where  the  scattered 
field  from  the  object  is  greater  than  the  numerical  er- 
rors introduced.  However , if  weak  scattering  events  are 
considered , such  as  for  objects  buried  in  the  ground , the 
dispersion  errors  can  be  of  the  same  order  of  magnitude 
as  the  scattered  field.  In  this  paper  we  derive  the  modi- 
fied Fresnel  coefficients  for  a homogeneous  lossy  ground 
which  are  consistent  with  the  FDTD  algorithm  and  can 
be  used  in  the  Huygens 7 sources  for  suppressing  numer- 
ical errors.  It  is  shown  that  the  modified  Fresnel  coeffi- 
cients reduces  the  noise  level  for  an  empty  FDTD  volume 
by  SO  to  60  dB.  As  an  application  example , the  far-zone 
scattering  of  a buried  dipole  is  considered. 

1 Introduction 

Numerical  simulation  of  electromagnetic  scattering  from 
objects  on  or  below  ground  is  an  important  topic  in  appli- 
cations such  as  ground  penetrating  radar  and  synthetic 
aperture  radar.  Both  time-domain  methods  and  fre- 
quency domain  methods  are  represented  among  recently 
published  work  [1-12].  Time-domain  methods  such  as 
FDTD  have  the  advantage  of  generating  results  for  the 
entire  frequency  range  in  a single  simulation.  In  scat- 
tering applications  where  the  transmitting  source  is  po- 
sitioned far  away  from  the  scatterer,  the  incident  field 
can  be  approximated  by  a plane  wave.  If  the  plane 
wave  is  created  using  a Huygens’  surface  enclosing  the 


object  [13],  a large  computational  volume  implies  large 
dispersion  errors  of  the  incident  field.  One  way  to  over- 
come this  for  free-space  simulations,  is  to  use  dispersion 
compensated  Huygens’  sources  [14],  applied  directly  in 
the  time  domain. 

When  a ground  is  included  as  a homogeneous  lossy 
half-space,  the  incident,  reflected  and  refracted  fields  are 
preferably  pre-calculated  in  the  frequency  domain  where 
it  is  relatively  simple  to  apply  the  Fresnel  reflection  and 
refraction  coefficients  to  the  incident  pulse  spectrum. 

However,  even  if  the  incident  fields  are  dispersion  com- 
pensated according  to  [14],  our  experience  shows  that 
the  errors  due  to  the  use  of  analytical  Fresnel  coefficients 
are  much  greater  than  the  dispersion  errors  occurring 
when  a plane  wave  propagates  across  the  computational 
volume.  The  reason  for  this  is  that  the  analytical  Fresnel 
coefficients  are  not  consistent  with  FDTD. 

Numerical  properties  of  reflection  coefficients  in  FDTD 
have  been  studied  in  several  papers.  Recently  Hirono 
et.  al.  [15]  derived  theoretical  reflection  coefficients  for 
a lossless  dielectric-dielectric  interface  for  a 2-D  case. 
In  [16]  reflection  coefficients  of  a staircased  air-dielectric 
interface  are  studied  for  a boundary  tilted  45°  with 
respect  to  the  FDTD-lattice.  Theoretical  expressions 
for  the  reflection  coefficients  at  interfaces  of  perfectly 
matched  layers  (PML)  have  also  been  presented  by  Fang 
et.  al.  [IT]. 

In  this  paper  we  derive  modified  FDTD  Fresnel  re- 
flection and  refraction  coefficients  for  oblique  incidence 
in  3-D  for  polarizations  parallel  as  well  as  perpendicular 
to  the  plane  of  incidence.  These  coefficients  are  intended 
for  plane  wave  excitation  in  scattering  applications  where 
the  ground  is  modeled  as  a homogeneous  lossy  dielectric 
half-space.  Using  these  Fresnel  coefficients,  significant 
improvements  of  scattering  results  for  weak  scatterers 
(such  as  buried  objects)  are  obtained. 

The  ground  model  used  here  is  a homogeneous  lossy 
dielectric  half-space  with  the  interface  aligned  with  the 
FDTD-lattice  axis.  The  tangential  electric  field  compo- 
nents in  the  top-layer  are  updated  using  the  constitutive 
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parameters  of  the  ground  while  the  normal  components, 
a half  cell  size  above,  are  entirely  updated  in  the  free- 
space  region.  Hence,  no  adjustment  of  the  interface  po- 
sition is  applied  using  harmonic  and  arithmetic  mean 
values  of  the  permittivities  at  the  interface  such  as  de- 
scribed in  [15]. 


2 Huygens’  sources  for  plane 
wave  incidence  on  a homoge- 
neous ground 


Plane  wave  incidence  upon  a homogeneous  half  space, 
using  a broad  band  pulse  requires  that  the  incident,  re- 
flected and  transmitted  fields  are  calculated  in  the  time- 
domain  on  a Huygens’  surface  enclosing  the  scattering 
object.  This  procedure  has  been  described  in  [1]  for 
stratified  media.  If  the  scatterer  is  large  it  is  normally 
favorable  to  use  the  total-field/scattered-field  method  us- 
ing a Huygens’  surface. 

The  principle  of  implementing  the  Huygens’  sources 
at  an  interface  between  free  space  and  a lossy  homoge- 
neous half-space  is  described  briefly  as  follows.  Above 
the  ground,  the  Fresnel  reflection  coefficients,  T y and 
r±,  are  calculated  in  the  frequency  domain  for  a given 
incidence  angle,  to  the  surface  normal.  The  incident 
pulse  spectrum  is  multiplied  with  the  corresponding  re- 
flection coefficient  and  the  result  is  transformed  into  the 
time-domain.  These  time-domain  fields  are  stored  in  ar- 
rays, one  for  each  field  component.  Given  the  retarded 
times  for  the  incident  and  reflected  fields  at  each  point  of 
the  Huygens’  surface,  the  fields  can  be  extracted  by  a ta- 
ble look-up  procedure  [18] ,[19].  The  situation  within  the 
homogeneous  ground  is  treated  in  a similar  way,  using 
the  Fresnel  refraction  coefficients  Ty  and  Tj_,  although  a 
lossy  ground  requires  storage  of  time-domain  field  com- 
ponents at  each  layer  in  the  ground  due  to  the  frequency 
dispersion. 

Consider  the  situation  in  Fig.  1.  The  analytical  Fres- 
nel coefficients  for  the  electric  field  perpendicular  and 
parallel  to  the  plane  of  incidence  are  given  by  [20] 
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Figure  1:  Reflection  and  refraction  at  an  interface  be- 
tween free  space  and  a homogeneous  half  space.  Left:  Po- 
larization perpendicular  to  the  plane  of  incidence.  Right: 
Polarization  parallel  to  the  plane  of  incidence. 


where 


is  the  index  of  refraction  of  the  lower  half-space,  9i  is 
the  incidence  angle  and  a non-magnetic  material  is  as- 
sumed. Note  that  for  a lossy  material  the  wave  is  an 
inhomogenous  plane  wave  and  the  complex  component 
Et,  for  parallel  incidence  in  Fig.  1,  has  no  real  direction. 
In  this  case  Ty  may  be  interpreted  as  the  relation  be- 
tween ( Zo/n)Ht  and  since  the  magnetic  fields  in 

both  media  are  parallel  (perpedicular  to  the  plane  of  in- 
cidence). This  will  be  used  in  section  3.2  when  discussing 
the  numerical  equivalence  to  Ty. 

Using  the  analytical  expressions  (1)  - (4)  in  FDTD 
gives  spurious  fields  at  the  Huygens’  surface  because  the 
FDTD-algorithm  is  not  consistent  with  the  analytical 
expressions.  Below  we  will  derive  the  Fresnel  coefficients 
that  are  consistent  with  the  FDTD-algorithm  and  we  will 
show  that  the  accuracy  can  be  increased  significantly  in 
scattering  simulations  of  weak  scatterers. 


3 Modified  Fresnel  coefficients 

To  derive  Fresnel  coefficients  consistent  with  FDTD,  the 
boundary  value  problem  at  the  interface  is  solved  by 
inserting  FDTD  consistent,  instead  of  analytical,  plane 
waves  into  the  FDTD  update  equations. 

The  derivations  below  are  divided  into  the  two  usual 
cases;  the  incident  electric  field  transverse  to  the  plane  of 
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incidence  (TE  or  _L)  and  the  incident  electric  field  paral- 
lel to  the  plane  of  incidence  (TM  or  ||).  A monochromatic 
plane  wave  of  amplitude  Uq  will  propagate  in  the  FDTD 
grid  according  to 


Due  to  (7),  the  direction  of  K — Kxx  + Kyy  + Kzz  will 
generally  not  coincide  with  k except  at  certain  angles. 
However,  the  difference  is  negligible,  see  Appendix  A. 
We  will  therefore  use 


U\lj)k  = U0etunAt-*r  (6) 

where  t/|£  ■ k is  an  electric  or  magnetic  field  component 
at  time  step  n at  position  (i,j>k)  in  the  FDTD  grid, 
k = kxx  + kyy  + kzz  is  the  numerical  wave  vector  and 
r = iAxx  + j Ayy  + kAzz. 

The  following  parameters  are  useful  for  the  deriva- 
tions. They  appear  naturally  when  substituting  the 
plane  wave  expressions  into  the  temporal  and  spatial  fi- 
nite differences  of  FDTD,  using  a time-average  for  the 
conduction  current  term; 
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where  ktz  is  the  z-component  of  the  numerical  wave  vec- 
tor in  the  ground,  which  is  complex  if  the  ground  is  lossy 
{a  7^  0).  Using  (7)  the  FDTD  dispersion  relation  in  free 
space  is  then  simply  written  as 


^ = kl  = kl  + kl  + kl  (8) 

co 

and  in  a homogeneous  lossy  ground  as 

Cr—  -jcr/Jow  = n2kl  = kl  + kl  + kfz  (9) 

co 

where  n is  the  FDTD  index  of  refraction 


Kx  = K0  cos  (f>  sin  6 (12) 

Ky  = Ko  sin  (j>  sin  0 
Kz  = KqcosO. 

which  will  actually  make  (11)  an  approximation.  Us- 
ing (9)  and  (12)  yields 

Ktz  = -K0 Vn2  - sin2  9 (13) 

where  we  have  chosen  the  minus  sign  since  it  is  assumed 
in  the  derivation  below  that  the  refracted  waves  will 
propagate  in  the  negative  ^-direction. 

Consider  the  situation  in  Fig.  1.  A plane  wave  is  in- 
cident on  a lossy  homogeneous  half-space  at  z < 0.  The 
plane  of  incidence  can  be  at  an  arbitrary  angle  0 with 
respect  to  the  horizontal  x-axis,  i.e.  we  do  not  restrict 
the  orientation  of  the  incident  plane  to  be  parallel  to  the 
x-  or  y-axes  in  the  FDTD-lattice.  Since  it  is  common  to 
express  the  Fresnel  coefficients  in  terms  of  the  incident 
angle  0*,  this  angle  will  be  used  below  instead  of  the  9- 
angle  for  the  incident  wave  vector.  Since  0i  = 7 r — 0, 
sin  0i  = sin  0 and  cos  9i  = — cos  0. 

3.1  Polarization  perpendicular  to  the 
plane  of  incidence  (TE) 

When  deriving  the  Fresnel  coefficients  for  the  TE-case, 
it  is  only  necessary  to  use  the  field  components  in  the 
X — Z- plane  shown  in  Fig.  2 (an  analysis  using  field 
components  in  the  Y — Z- plane  yields  the  same  results). 
The  exact  position  of  the  interface  between  the  upper 
and  lower  half-spaces  - illustrated  with  a faded  shading 
in  Fig.  2 - is  so  far  undetermined.  The  only  assumptions 
are  that  the  Hx\ij ^+1/2  and  ^ykj.fc+i  components  in 
Fig.  2 are  in  the  upper  half-space  and  that  the  others 
are  in  the  lower  half-space.  This  is  indicated  with  the 
horizontal  dashed  line  in  Fig.  2. 

We  also  need  to  define  the  different  horizontal  electric 
field  components  Ex  and  Ey  for  the  incident,  reflected 
and  refracted  waves.  In  the  upper  medium  we  have 


~2  • & 
n = er  - j — r- 

60W 

(10) 
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— sin  (j)€?ti)t-3k*z-3kyy-3ktz 

(14a) 

The  incident  wave  vector  can  be  written  as 

II 
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The  update  equation  for  Ey\ijtk  within  the  ground  is 


€r€0  ~ | -f  cr  1 


,n  — 1/2  ' 


in+l/2  n in— 1/2  ' 


At 

/ 

A z 
Ax 


(19) 


Inserting  the  plane  wave  expression  (17)  for  Ey  into  (19) 
the  left  hand  side  (LHS)  of  (19)  becomes 


LHS  = (ere0 joj  + a)Tx  cos  <j).  (20) 


Figure  2:  FDTD  field  components  used  in  the  derivation 
of  Fresnel  coefficients  for  the  TE-case.  The  dashed  line 
does  not  represent  the  interface  position , it  just  divides 
the  field  components  into  the  two  different  half-spaces . 


The  magnetic  field  component  Hx\ij,k+i/2  on  the  right 
hand  side  of  (19)  is  in  the  free  space  region  and  there- 
fore consists  of  an  incident  and  a reflected  part.  Us- 
ing (15),  (18a)  and  (18b)  we  get 

Hx\i,j,k+i/2  = - ( ' e -jil££  - fxe7^)  cos$!>  (21) 


and  in  the  lower  medium  we  have 


The  other  Hx  component  in  (19)  is  in  the  lower  half- 
space, using  (17)  and  (18c)  yield 


El  = —T±  sin  (16) 

El  = fxcos  (17) 

where  T±  and  T±  are  the  FDTD  Fresnel  coefficients  yet 
to  be  determined.  We  also  need  the  relationship  between 
the  electric  and  magnetic  fields  for  each  of  the  plane 
waves.  Utilizing  (7)  we  can  rewrite  the  FDTD  update 
equations  for  some  of  the  magnetic  fields  as 


HojuHlx 
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Hx\i,jik-1/2  = ~T±e^  cos0.  (22) 

PqOJ 

The  Hz -components  are  in  the  lower  medium.  Us- 
ing (16),  (17)  and  (18d),  the  spatial  difference  of  Hz 
in  (19)  can  then  be  written  as 


Hz  \i+l/2,j,k  Hz  [i— 1/2, j,/: 
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~ K ~ ~ 

= -jT± — %:  (Kx  cos  + Ky  sin  </>) 
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where  we  also  have  used  (12).  Using  (20),  (21),  (22) 
and  (23)  in  (19)  now  yield 


By  inserting  the  plane  wave  expressions  (14)  - (17) 
into  the  FDTD  update  equations  for  Ey\ijtk  and 
H*  li,  jtk+ 1/2 j see  Fig.  2,  we  are  able  to  determine  the 
unknown  coefficients  T±  and  T±.  We  choose  the  phase 
reference  point  at  the  position  of  Ey\ij7k  in  the  middle 
of  Fig.  2.  Note  that  all  components  in  Fig.  2 have  the 
same  ^/-coordinate  (y  = 0). 
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The  other  equation  needed  is  the  update  equation  for 
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above  the  ground: 
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Inserting  plane  wave  expression  for  the  Hx-  and  Ey - 
fields  using  (15),  (17),  (18a)  and  (18b)  yields 

-jKz  (e~jL^  - f ±ejL£i')  = 

( e-jk;Az  + fxej~k*Az  \ (26) 

( j ~ ~K~z 

Solving  the  equation  system  (24)  and  (26)  with  respect 
to  the  two  unknowns  T j_  and  T±  yields 


a cos  0{  — fi\jh2  — sin2  6i 

a*  cos  6i  + /3  \Jn2  — sin2 
(a  4-  cm*)  cos#* 
a*  cos  Oi  -F  fi  \/n2  — sin2  #* 


(27)  Figure  3:  FDTD  field  components  used  in  the  derivation 
of  Fresnel  coefficients  for  the  TM-case.  The  dashed  line 

(28)  does  not  represent  the  interface  position , it  just  divides 
the  field  components  into  the  two  different  half-spaces. 


where  we  have  used  (7c),  (7d),  (12)  and  (13).  The  coeffi- 
cients q and  fi  are  factors  corresponding  to  propagation 
across  a half  cell  in  the  ^-direction  in  the  upper  and 
lower  medium  respectively 
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Note  that  the  only  ^-dependence  of  (27)  and  (28)  are 
through  the  FDTD  dispersion  relation  and  the  numerical 
wavenumber  in  (29)  and  (30). 

As  expected,  (27)  and  (28)  reduces  to  the  analytical 
equations  (2)  and  (4)  as  Az,  Ay,  Az  — > 0.  Also,  note  that 
when  |n|  ->  oo,  as  for  a perfect  conductor,  T±  — 1, 
T]_  — > 0,  i.e.  the  reflection  occurs  at  the  layer  of  electric 
field  components  at  z = 0,  while  for  n = 1 (when  a = fi) 

r±  = o,  fx  = i.  _ 

If  ko\h\Az  C 1,  rx  « a2rx  and  Tx  as  ot/fiT±,  which 
are  equal  to  the  analytical  coefficients  with  the  interface 
at  z — Az/2,  as  shown  in  Appendix  B. 


3.2  Polarization  parallel  to  the  plane  of 
incidence  (TM) 

The  coefficients  for  the  TM-case  can  be  found  in  a similar 
manner.  Consider  the  components  in  the  X — Z -plane 
in  Fig.  3.  In  this  case  only  the  Hz -components  are  zero. 
As  in  the  TE-case  we  only  need  to  consider  the  field 
components  at  y = 0.  In  the  upper  medium  we  have 
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(33a) 

Erz  = f||  sin0iejut-j^x~j'kyy+j'klZ 
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and  in  the  lower  medium  we  have 

El  = 7^  cos 

(34) 

Ely  = f||h  sin 

(35) 

El  = T^ve3ut-jkxx-jkyy-jkttZ 

(36) 

where  Ty  is  the  reflection  coefficient  and  and  Tyv 
are  horizontal  and  vertical  refraction  coefficients,  respec- 
tively. 
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We  also  need  relationships  between  Ex  and  Hy  for  the 
plane  waves  in  both  media  analogous  to  (18): 


Solving  (39)  and  (41),  using  (12)  and  (13),  with  respect 
to  the  unknowns  Fy  and  Tp  yields 


ju>e0E'x  = 

(37a) 

~ /? n2  cos  0i  — ay/ n2  — sin2  0i 

jue0Erx  = 

-jKzHy 

(37b) 

!>  fin2  cos  0i  4-  a*  y/n2  — sin2  0i 

(juie0er  + a)Etx  = 

jKtzHy 

(37c) 

~ (a  4-  a*)  cos#*\/n2  — sin2  0i 

1 ~ — 

The  update  equation  for  Hy |, 

,j,k+ 1/2  in  Fig.  3 is 

0n2  cos  0i  4-  a*  y/n2  — sin2  0i 

-Mo 


' it  in+1  _ ij  In 

ny\i,j,k+if2 

A t 


/ VP  I w+l/2  rp  in+1/2 

Az 


Finally  we  need  to  determine  T\\v.  Using  V • E£ 
within  the  lower  medium 

Kt  E*  = T]{hK0  sinOi  - T\\vK0\ff>?AtiP  0t  = 0 
(38)  which  together  with  (43)  yields 


(42) 

(43) 
= 0 

(44) 


/ rp  in+1/2  _ T?  |n_{_1/2 


T\u,  = 


(a  4-  a*)  sin#*  cos#* 

/?n2  cos  4*  a*  y/h2  — sin2  #* 


(45) 


We  now  define  the  phase  reference  point  at  the  position 
of  Ex\ijfk  (i*e.  the  same  z-level  as  for  TE-polarization), 
the  only  component  in  (38)  that  is  in  the  lower  medium, 
see  Fig.  3.  Inserting  the  plane  wave  expressions  (31)  - 
(36)  into  (38)  and  using  (37)  gives 


As  in  the  analytical  case  we  may  define  a total  T\\  as  the 
relation  between  (Z0/n)H£  and  ZqH1  (which  are  paral- 
lel), where  H = (Kx  E)/d;/io-  This  yields 

- _ («  + <«>  cos  a,  (46) 

fin2  cos  9i  + a*  \/n2  — sin2  Oi 


UJ€q  .:  k i A 2 ./  kji  Ai  \ n , 

-ju>no[-^e  1 2 H — =^-rneJ  2 )cos0jcos0 
. Kz  Kz 

' e-jk,Az  _f  ejk,A 


cos  Qi  cos  <j> 


Tjjftcos^ 


Az 


— / • ; k z A z \ 

4-  jKx  sin  f e 72  4-Fn^  2 J 


(39) 


The  other  update  equation  needed  is  for  Ex\ij^  in  the 
lower  medium. 


En 


CrCo 


in+1/2 


■E2 


in— 1/2 ' 
\iijyk 


At 


+ cr 


in+1/2 

Uj>k 


+ EX  | 
2 


n— 1/2 


When  Ax,Ay,Az  — ¥ 0,  the  numerical  Fresnel  coeffi- 
cients (42)  and  (46)  converge  to  (1)  and  (3),  respectively. 
Also,  as  expected,  when  Oi  = 0,  Tp  = 0,  Tp  — T±  and 
r||  = — T+.  The  behaviour  when  ko\h\Az  4C  1 is  similar 
to  the  TE-case  although  the  phase  variation  is  slightly 
more  complicated,  see  Appendix  B. 

The  derived  Fresnel  coefficients  contain  the  FDTD  in- 
dex of  refraction  (10)  consistent  with  the  FDTD  disper- 
sion relation  (9),  and  they  are  also  modified  with  the  fac- 
tors a and  ft  compared  to  the  analytical  ones.  They  are 
therefore  referred  to  as  “modified”  Fresnel  coefficients  in 
the  following  text. 


' u in  _ u in 

nV\iyj,k+l/2 

Az 


1/2 


(40) 


Inserting  the  plane  wave  expression  (31)  - (36)  into  (40) 
and  using  (37)  gives 


(ere0 + <?)T\\h  cos  (j>  = - e cos  0*  cos  <j> 

WeoT||  j jjf,;,  Ar  n I 

— cosq. cos© 

KzAz 

(Cr^Ojd)  4-  v)T\\h  • Z&.L  , 

H sr -~eJ^  cos  6 

jKtzAz 


(41) 


4 Reflection  and  refraction  coeffi- 
cients in  a 1-D  FDTD  case 

To  demonstrate  the  accuracy  of  the  model,  a one  dimen- 
sional FDTD  simulation  was  performed  where  a Gaus- 
sian pulse  was  incident  on  a homogeneous  medium  with 
er  = 10  and  a = 0.01  S/m.  The  reflected  tangential  77- 
field  and  the  refracted  tangential  E-field  adjacent  to  the 
interface  were  registered.  The  results  were  transformed 
into  the  frequency  domain  and  the  reflection  and  refrac- 
tion coefficients  were  calculated.  The  results  are  shown 
in  Fig.  4 together  with  the  analytical  Fresnel  coefficients. 
The  cell  size  was  0.01  m and  the  maximum  frequency  2 
GHz  corresponds  to  15  cells /wavelength  in  free  space  and 
only  five  cells/ wavelength  in  the  lossy  medium. 
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Figure  4:  Reflection  and  refraction  at  normal  incidence 
on  an  interface  between  free  space  and  a homogeneous 
half  space.  Left:  Reflection  coefficient.  Right:  Refrac- 
tion coefficient.  Note  that  the  FDTD  results  and  the 
predictions  according  to  (27)  and  (28)  are  iu distinguish- 
able. 


Figure  5:  X-Z  slice  of  the  computational  volume  at  y — 
20 A y.  The  Huygens'  surface  extends  10  cells  below  the 
ground.  The  far-zone  transformation  is  performed  in  the 
scattered  field  region  2 cells  outside  the  Huygens'  surface. 


As  seen  from  Fig.  4,  the  correspondence  between  the 
predicted  FDTD  results  and  the  simulated  FDTD  results 
are  excellent. 

5 Plane  wave  incidence  on  an 
empty  3-D  FDTD  volume 

The  modified  Fresnel  coefficients  (27),  (28),  (42),  (43) 
and  (45),  were  implemented  in  a 3-D  FDTD  code  where 
they  were  used  in  the  Huygens’  routine  as  well  as  in  the 
near-  to  far-zone  transformation  routine  [12].  Scattering 
simulations  were  performed  for  an  empty  FDTD-volume 
in  order  to  establish  the  lower  limit  of  Radar  Cross  Sec- 
tion (RCS)  simulations  for  a given  problem  size.  Simu- 
lations were  also  performed  using  the  analytical  Fresnel 
coefficients  as  a reference  case.  The  plane  of  reflection 
for  the  analytical  coefficients  (l)-(4)  was  positioned  at 
the  layer  of  tangetial  magnetic  field  components,  Az/2 
above  the  phase  reference  point  used  in  the  derivations 
of  the  modiefied  Fresnel  coefficients.  This  gives  the  most 
accurate  results  for  analytical  Fresnel  coefficients.  A mo- 
tivation for  this  can  be  found  in  Appendix  B. 

When  calculating  the  Huygens’  sources  in  the  free- 
space  region,  the  incident  and  reflected  pulse  propaga- 
tions were  dispersion  compensated  according  to  [14].  For 
the  Huygens’  sources  within  the  ground,  this  dispersion 
compensation  procedure  was  applied  on  the  propagation 


occurring  in  the  free  space  region  before  entering  the 
lossy  ground.  The  dispersion  of  the  pulse  due  to  propa- 
gation within  the  lossy  ground  was  then  taken  into  ac- 
count, using  the  numerical  wavenumbers  when  calculat- 
ing the  pulses  in  the  frequency  domain  at  each  horizontal 
layer.  The  dispersion  compensation  was  not  used  in  the 
reference  case  since  its  effect  is  small  compared  to  the 
adjustment  of  the  Fresnel  coefficients. 

The  pulse  functions  were  pre-calculated  and  stored 
in  arrays,  field  values  were  extracted  using  a table-look 
up  procedure  [18]  and  a third  order  Lagrange  interpo- 
lation. A Gaussian  derivative  pulse  was  used  for  the 
incident  plane  wave.  The  pulse  function  was  Einc  = 
ry/2eexp(— r2)  V/m,  where  r = 4/(/3At)(t  — 1.5/3 At) 
with  /?  = 80.  The  computational  volume  consisted  of 
50  x 40  x 50  cubic  cells,  see  Fig.  5 ( Ax  = Ay  = Az  = 0.01 
m).  The  outer  boundary  was  terminated  by  a 6-cell-thick 
PML  [21].  The  PML  within  the  ground  was  treated  in  a 
similar  way  as  in  [22]. 

A Huygens’  surface  was  positioned  five  cells  from  the 
outer  boundary.  The  upper  level  of  the  lossy  half-space 
was  15  cells  from  the  lower  PML-region,  see  Fig.  5.  As 
before,  the  parameters  for  the  ground  were  er  = 10  and 
a = 0.01  S/m  (A  = 15  cells  at  2 GHz  in  the  free  space 
region  and  only  5 cells  within  ground). 

Three  different  incident  angles  were  used  in  the  sim- 
ulations; 6i  — 0° (normal  incidence),  0*  = 45°,  (j)  = 0 
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Figure  6:  Monostatic  RCS  of  empty  FDTD  volume  as 
function  of  frequency.  Incidence  angle  $i  = 0°. 


and  = 45°,  (p  — 30°.  Both  TE-  and  TM-  polariza- 
tions were  used  in  the  = 45°,  <f>  = 0°-case.  In  the 
other  cases  TM-polarization  were  used  (both  polariza- 
tions coincide  for  the  0i  — 0°-case).  The  scattered  field 
was  determined  using  a time-domain  near-zone  to  far- 
zone  transformation  applied  three  cells  from  the  outer 
boundary  [12].  Both  monostatic  and  bistatic  RCS  results 
were  calculated.  The  results  using  the  modified  Fresnel 
coefficients  are  labeled  “Modified”  in  Figs.  6-  12.  The 
corresponding  results  using  the  analytical  expressions  for 
the  Fresnel  coefficients  (including  phase  corrections)  are 
labeled  “Reference”. 

The  effects  of  using  modified  Fresnel  coefficients  on  the 
far-zone  scattered  field  are  illustrated  in  Figs.  6-9  where 
the  monostatic  and  bistatic  RCS  of  an  empty  FDTD  vol- 
ume are  shown  for  four  cases  of  incident  directions.  The 
scattered  field  is  calculated  in  the  same  polarization  as 
the  incident  field.  The  volume  enclosed  by  the  Huygens’ 
surface  is  0.4  x 0.3  x 0.4  m3. 

As  seen  from  the  figures,  the  RCS  level  of  the  empty 
volume  has  been  reduced  by  30  - 60  dB  using  the  mod- 
ified Fresnel  coefficients.  The  RCS  level  is  higher  for 
TM-polarization  than  for  TE-polarization  in  the  refer- 
ence case,  according  to  Figs.  7 and  8.  The  reason  is  that 
the  reflection  plane  offset  of  A^/2  is  a better  approxi- 
mation for  TE-  than  for  TM-polarization  according  to 
Appendix  B. 

The  use  of  modified  Fresnel  coefficients  does  not  af- 
fect the  computational  time  during  the  FDTD  main  loop 
since  the  incident  pulse  functions  are  pre-calculated  and 
stored  in  tables.  The  time  increase  in  using  the  modified 
Fresnel  coefficients  compared  to  the  analytical  ones  has 
been  measured  to  « 5%  when  calculating  the  frequency 


Figure  7:  Monostatic  RCS  of  empty  FDTD  volume  as 
function  of  frequency.  Incidence  angle  $i  = 45°,  0 = 0°. 
Both  TE-  and  TM - results  are  shown . 


domain  pulse  spectra  before  the  main  loop.  In  the  near- 
to  far-zone  transformation,  the  Fresnel  coefficients  are 
applied  after  the  main  loop  according  to  [12]  which  con- 
sequently does  not  affect  the  computational  time  of  the 
main  loop  either. 

The  demonstrated  reduction  of  the  RCS  noise  level  is 
important  when  studying  low  level  scattering,  such  as  for 
buried  objects.  An  example  of  this  follows  in  section  6. 

6 Scattering  from  a dipole 

A resonant  dipole  buried  0.2  meters  below  the  ground 
level  was  chosen  as  an  example  of  a scattering  object. 
The  dipole  length  was  0.21  m and  the  radius  1 mm.  The 
RCS  of  the  dipole  was  calculated  using  the  Method  of 
Moment  program  NEC-3  and  the  FDTD  method,  both 
with  analytical  and  modified  Fresnel  coefficients.  A thin 
wire  model  according  to  [18]  was  used  for  the  dipole 
which  was  oriented  in  the  ^-direction.  The  dipole  was 
modeled  using  20  cells,  which  corresponds  to  a true 
dipole  length  of  0.21  m since  the  staircasing  effect  in 
FDTD  causes  the  electric  dipole  length  to  increase  with 
a half  cell-size  at  each  end. 

The  computational  volume,  cell  size,  positions  of  Huy- 
gens’ and  near-  to  far-zone  transformation  surfaces  and 
ground  parameters  were  identical  as  in  the  previous  sec- 
tion. The  only  difference  was  the  position  of  the  ground 
level,  which  was  35  cells  from  the  lower  PML-boundary, 
see  Fig.  10.  The  dipole  was  positioned  0.205  m below  the 
ground  level  in  the  NEC-3  simulation  since  the  main  re- 
flection occurs  at  the  level  of  tangential  magnetic  fields, 
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Figure  8:  Bistatic  RCS  of  empty  FDTD  volume  as  func- 
tion of  frequency.  Incidence  angle  6{  = 45°,  0 = 0°. 
Scattering  in  specular  direction.  Both  TE-  and  TM-  re- 
sults are  shown . 


as  discussed  in  Appendix  B. 

The  plane  wave  was  incident  in  the  same  plane  as  the 
dipole  at  0*  = 45°.  The  frequency  band  was  restricted 
to  be  below  1 GHz  using  a slightly  narrower  pulse  spec- 
trum, compared  to  the  empty  volume  simulations  in  sec- 
tion 5 (er  = 10  and  a — 0.01  S/m  yield  approximately 
10  cells/ wavelength  at  1 GHz  in  the  ground).  The  inci- 
dent field  was  applied  using  a modulated  Gaussian  pulse 
with  Einc  = cos(27r/o[t  ~ 1.5/?A£])  exp(— r2)  V/m,  where 
r = 4/(/?A t)(t  — 1.5/3  At)  with  ft  = 320  and  /0  = 0.5 
GHz. 

As  a reference,  the  monostatic  RCS  was  calculated 
without  the  dipole  in  FDTD  using  both  the  modified  and 
analytical  Fresnel  coefficients  in  the  Huygens  routine. 

The  monostatic  RCS  of  the  dipole  can  be  seen  in 
Fig.  11.  Both  the  modified  and  the  analytical  Fresnel  co- 
efficients have  been  used  in  the  FDTD  simulations.  For 
comparison,  the  scattering  result  for  an  empty  FDTD 
-volume  without  the  dipole  using  the  analytical  Fresnel 
coefficients  is  also  shown  in  the  graph  as  a dashed-dotted 
line  (the  RCS-level  of  the  empty  FDTD-volume  using 
the  modified  Fresnel  coefficients  are  below  -120  dBm2 
in  this  case).  As  before,  results  using  analytical  Fres- 
nel coefficients  are  labeled  “Reference”.  Obviously,  the 
scattered  field  in  the  reference  case  without  the  dipole  is 
of  the  same  order  as  the  scattered  field  from  the  dipole. 
Hence,  the  poor  correspondence  with  the  NEC-result  in 
this  case.  There  is  a difference  between  the  NEC-  results 
and  the  FDTD-results  using  modified  Fresnel  coefficients 
which  is  due  to  the  poor  resolution  in  the  ground  (only 
10  cells  per  wavelength  at  1 GHz)  and  to  the  thin-wire 


Figure  9:  Monostatic  RCS  of  empty  FDTD  volume  as 
function  of  frequency.  Incidence  angle  0j  = 45°,  (j)  = 30°. 
TM-polarization . 


model  in  FDTD.  Also,  since  the  modified  Fresnel  coeffi- 
cients are  different  from  the  ones  in  NEC-3,  the  refracted 
fields  reaching  the  dipole  are  slightly  different  in  the  two 
models.  This  has  been  verified  in  a simulation  using 
twice  as  high  resolution  in  FDTD,  which  showed  better 
convergence  between  FDTD-  and  NEC-  results. 

It  is  often  of  interest  to  study  the  time-domain  re- 
sponse of  the  electric  fields  for  these  type  of  scatter- 
ing simulations.  The  scattered  electric  far-fields  from 
the  NEC-3  were  therefore  Fourier  transformed  into  the 
time-domain  after  multiplication  by  the  incident  pulse 
spectrum  used  in  the  FDTD-simulations.  The  far-fields 
in  FDTD  were  calculated  in  the  time-domain  using  the 
method  described  in  [12].  The  phase  reference  point  was 
set  at  the  surface  straight  above  the  dipole,  see  point 
P in  Fig.  10.  To  get  a distance  independent  result,  the 
voltage  rEe  was  extracted  from  the  simulations. 

FDTD  results  without  the  dipole  is  shown  in  top  of 
Fig.  12.  using  both  the  modified  and  analytical  Fres- 
nel coefficients.  In  the  reference  case,  the  signal  appears 
slightly  before  t = 0 due  to  the  disturbance  at  left  inter- 
section between  the  ground  and  the  Huygens1  surface  in 
Fig.  10. 

The  time  domain  results  with  the  dipole  in  place  is 
shown  in  the  bottom  of  Fig.  12  including  the  NEC-3 
result.  As  seen  from  the  graphs,  the  correspondence  is 
very  good  except  for  the  FDTD  reference- case.  After  the 
disturbance  from  the  Huygens’  surface  has  vanished,  the 
signals  are  nearly  identical.  If  a larger  region  would  be 
used,  the  error  would  increase  as  the  side  length  of  the 
Huygens’  surface  increases. 
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Figure  10:  X-Z  slice  of  the  computational  volume  at  y = 
20 Ay  for  the  dipole  case.  The  Huygens'  surface  extends 
30  cells  below  the  ground.  The  far-zone  transformation  is 
performed  in  the  scattered  field  region  2 cells  outside  the 
Huygens'  surface.  The  dipole  is  positioned  20  cells  below 
the  ground  level  (of  tangential  electric  field  components 
in  FDTD) 


7 Conclusions 

Fresnel  coefficients  for  a homogeneous  lossy  dielectric 
half-space,  consistent  with  the  FDTD  algorithm  have 
been  derived.  The  coefficients  are  valid  for  arbitrary  inci- 
dence angles  in  a 3D  FDTD  simulation.  Using  these  co- 
efficients in  the  Huygens’  routine  reduce  the  errors  signif- 
icantly. The  improvements  achieved  have  been  demon- 
strated for  a scattering  simulation  of  a buried  dipole  in 
a homogeneous  lossy  dielectric  half-space.  The  modified 
Fresnel  coefficients  lowers  the  limits  of  RCS-levels  possi- 
ble to  simulate  with  FDTD,  where  plane  wave  incidence 
on  a homogeneous  ground  is  considered. 


Figure  11:  Monostatic  RCS  of  buried  dipole  0.205  m 
below  the  ground  level.  Incidence  angle  0i  = 45°,  </>  — 0°. 
TM-  polarization. 


Time  (ns) 


Figure  12:  Time  response  of  the  monostatic  scattering 
from  a buried  dipole  0.205  m below  the  ground  level  (bot- 
tom) and  without  dipole  (top).  Incidence  angle  0*  — 45°, 
0 — 0°.  TM-polarization. 
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Appendix  A Incident  angle  ap- 
proximation 


The  validity  of  the  numerical  incident  angle  approxima- 
tion can  be  illustrated  by  taking  the  dot  product  between 

k = k/k0  and  K = K/AV  Using  the  series  expansion 
sin  a;  « x — x3/6  + 0(x 5)  in  (7a)-(7c),  the  dot  product 
can  be  written  as 


~ * k • K kr.Kr  4*  kyKv  + k+Kz 

k ■ K = 11  - — — — 


ko  A"o 

jfc2  _ j_  i.2 


24  + 


k0I<o 

K*  y7 


koy/K*  + K*  + K* 
+ 0 [Ckv^r}  ,u  = x,y,z. 


_____  _j_  JU2  __ 

24  24  (A.l) 


The  square  root  can  be  expanded  further  and  if  we  in- 
troduce the  variable  p,  where 


p2  = s4xAx2  + s^Ay2  + siAz2  (A.2) 


= cos  <^sin  9,  sy  — sin  0sin  0,  sz  ~ cosO  (A. 3) 


we  can  rewrite  (A.l)  as 


k-K 


H [^1  - ^§4")  + 0(fcoP4)] 
^ + 0(^oP4) 


(i  - + 0(Hp4) 

(i  - %-)  + 0(k$p*) 


1 + O(kop4). 


(A.4) 


This  result  indicates  that  the  angle  approximation  (12) 
is  very  accurate. 


where  we  have  introduced  N = y/n2  — sin2  0*  and  used 

ktz  ~ —koN 

kz  — —ko  cos  0i.  (B.2) 

Rewriting  (B.l)  yields 

r(cos^-jv)  + 

(cos  0|  + at)  + 

= a2r±  + O [(fcoiVAz)2]  . 

(B.3) 

This  indicates  that  for  low  frequencies  -and  not  too  high 
|n|-  the  modified  reflection  coefficient  Fj_  converges  into 
the  analytical  reflection  coefficient  times  the  phase  factor 
a2,  corresponding  to  a shift  in  the  reference  point  by 
Az/2  in  the  ^-direction.  A similar  analysis  for  gives 

Tjl  = + O [(fcoiVAz)2]  (B.4) 

which  also  corresponds  to  a shift  of  the  reference  point 
by  Az/2  (a)  and  an  attenuation  of  the  plane  wave  for 
that  distance  (1//?). 

Analogous  analysis  for  T |j  and  Tjj  gives  slightly  more 
complicated  results.  It  can  be  shown  that, 

f ||  = a2+-y  r,|  + o [(fcoiVAz)2]  (B.5) 

and 

r^^TH+O^oiVA,)2]  (B.6) 

where  7,  eq  and  62  are  proportional  to  for  small  values 
of  6i.  These  phase  factors  are  thus  not  so  well  approx- 
imated by  a simple  shift  of  the  reference  plane  as  for 
perpendicular  (TE)  polarization. 


+ O [(fcoNAz)2] 
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